% function [FEAS_FLAG, bu_a, info] = minPICU_mt(H, bu_min, w, cb, opts)
% minPICU_mt
%
% Multitone extension of minPICU for a scalar U-user Gaussian IC.
%
% INPUTS
%   H       U x U x N channel, where H(r,u,n) is user u to receiver r on tone n.
%           A U x U matrix is accepted and treated as N=1.
%   bu_min  U x 1 requested total user rates, summed over tones.
%   w       U x 1 positive user energy weights.
%   cb      1 for complex baseband, 2 for real baseband.
%   opts    same style as minPICU; opts.useParfor optionally parallelizes tones.
%
% ATOMS
%   Uses the IC atomization:
%       atoms (u,S) with u in S.
%       K = U * 2^(U-1).
%   For U=2 this gives 4 atoms.
%   For multitone operation, powers and rates are K x N arrays.
%
% IMPORTANT FIRST-VERSION DESIGN CHOICE
%   This version uses one receiver order set per harvested vertex and applies
%   that order set on every tone. The power/rate optimization is tone-separable
%   for a fixed theta and order set, and the outer theta is common across tones.
%   This is the safest first multitone step because it preserves the scalar
%   minPICU machinery while adding per-tone loading and parallelism.
function [FEAS_FLAG, bu_a, info] = minPICU_mt(H, bu_min, w, cb, opts)
    if nargin < 5
        opts = struct();
    end
    opts = mp_set_defaults(opts);

    if isfield(opts, 'randomSeed')
        rng(opts.randomSeed);
    end

    if ndims(H) == 2
        H = reshape(H, size(H,1), size(H,2), 1);
    end

    [Ur, Uc, N] = size(H);
    assert(Ur == Uc, 'H must be U x U x N.');
    U = Ur;

    bu_min = bu_min(:);
    w = w(:);

    assert(numel(bu_min) == U, 'bu_min must have length U.');
    assert(numel(w) == U && all(w > 0), 'w must have length U and be positive.');
    assert(cb == 1 || cb == 2, 'cb must be 1 or 2.');

    G = abs(H).^2;
    atoms = mp_build_atom_table(U);
    K = numel(atoms);

    if opts.verbose
        fprintf('minPICU_mt:\n');
        fprintf('  U = %d, K = %d atoms, N = %d tones\n', U, K, N);
        fprintf('  target = [%s]\n', num2str(bu_min.', ' %.6g'));
        fprintf('  weights = [%s]\n', num2str(w.', ' %.6g'));
        fprintf('  cb = %d\n', cb);
        fprintf('  IC atomization (own receiver always included), exact Imin reward\n');   
    end

    if ~isempty(opts.thetaMaxOverride)
    theta_max = opts.thetaMaxOverride(:);
      else
      theta_max = mp_compute_theta_max_bounds(G, bu_min, w, cb, opts, atoms);
    end
    theta_max = max(theta_max(:), opts.thetaFloor);
    theta0 = 0.5 * theta_max(:);

    xc = theta0(:);
    A = diag(theta_max(:).^2);

    if opts.verbose
        fprintf('  theta_max = [%s]\n', num2str(theta_max.', ' %.6g'));
        fprintf('  theta0    = [%s]\n', num2str(theta0.', ' %.6g'));
    end

    theta_history = zeros(U, opts.maxEllipIter);
    ellip_trace = repmat(struct( ...
        'iter', [], ...
        'theta', [], ...
        'nOrders', [], ...
        'best_bu', [], ...
        'best_bu_support', [], ...
        'best_E', [], ...
        'shortfall', [], ...
        'pool_idx', []), opts.maxEllipIter, 1);

    vertex_pool = mp_init_vertex_pool(K*N, U);
    prev_best = [];

    for it = 1:opts.maxEllipIter
        theta_history(:, it) = xc;

        orderSet = mp_generate_order_candidates(xc, atoms, U, opts, prev_best);

        bestIter = [];
        idx_kept = [];

        for kk = 1:numel(orderSet)
            cand = mp_inner_solve_for_orders(G, bu_min, w, cb, opts, atoms, orderSet{kk}, xc, prev_best);
            [vertex_pool, idxAdded] = mp_pool_add_vertex(vertex_pool, cand, opts.poolUniqueTol);
            idx_kept = [idx_kept; idxAdded]; %#ok<AGROW>

            if isempty(bestIter) || cand.lag_obj < bestIter.lag_obj
                bestIter = cand;
            end
        end

        prev_best = bestIter;
        shortfall = bu_min(:) - bestIter.bu(:);

        if opts.thetaAdaptiveExpand
            hit = (xc >= opts.thetaNearMaxFrac * theta_max) & ...
                  (shortfall > opts.thetaExpandShortTol);

            if any(hit)
                old_theta_max = theta_max;

                % IMPORTANT: adaptive "expansion" must be monotone.
                % The earlier version used min(opts.thetaMaxCap, factor*theta_max),
                % which SHRANK entries when thetaMaxCap was smaller than the
                % current bound, e.g. 4.8e7 -> 1e6.  That invalidates the
                % ellipsoid box and can starve the affected users.
                theta_candidate = theta_max;
                theta_candidate(hit) = opts.thetaExpandFactor * theta_max(hit);

                if isfinite(opts.thetaMaxCap)
                    theta_candidate(hit) = min(opts.thetaMaxCap, theta_candidate(hit));
                end

                % Enforce no decrease, even if the user supplied a too-small cap.
                theta_max = max(old_theta_max, theta_candidate);

                if any(theta_max < old_theta_max - 100*eps(max(1,norm(old_theta_max,inf))))
                    error('theta_max decreased: old=[%s], new=[%s]', ...
                        num2str(old_theta_max.', ' %.6g'), ...
                        num2str(theta_max.',     ' %.6g'));
                end

                for uu = 1:numel(theta_max)
                    if theta_max(uu) > old_theta_max(uu)
                        A(uu,uu) = max(A(uu,uu), theta_max(uu)^2);
                    end
                    xc(uu) = min(xc(uu), theta_max(uu));
                end

                if opts.verbose
                    if any(theta_max > old_theta_max)
                        fprintf('    expanded theta_max -> [%s]\n', ...
                            num2str(theta_max.', ' %.6g'));
                    else
                        fprintf('    theta_max expansion requested but already at/beyond cap; unchanged -> [%s]\n', ...
                            num2str(theta_max.', ' %.6g'));
                    end
                end
            end
        end

        ellip_trace(it).iter = it;
        ellip_trace(it).theta = xc(:);
        ellip_trace(it).nOrders = numel(orderSet);
        ellip_trace(it).best_bu = bestIter.bu(:);
        ellip_trace(it).best_bu_support = bestIter.bu_support(:);
        ellip_trace(it).best_E = bestIter.E;
        ellip_trace(it).shortfall = shortfall(:);
        ellip_trace(it).pool_idx = idx_kept(:).';

        if opts.verbose
            fprintf(['  iter %2d: theta=[%s], orders=%d, E=%-.6g, ' ...
                     'bu=[%s], bu_sup=[%s], short=[%s]\n'], ...
                it, num2str(xc.', ' %.5g'), numel(orderSet), bestIter.E, ...
                num2str(bestIter.bu.', ' %.5g'), ...
                num2str(bestIter.bu_support.', ' %.5g'), ...
                num2str(shortfall.', ' %.3e'));
        end

        if norm(shortfall, inf) <= opts.hitTol
            break;
        end

        g = shortfall(:);
        if norm(g) <= 1e-14
            break;
        end

        denom = sqrt(g.' * A * g);
        if denom <= 1e-14
            break;
        end

        ntheta = U;
        Ag = A * g;
        xc = xc + (1/(ntheta+1)) * (Ag / denom);
        A = (ntheta^2/(ntheta^2-1)) * ...
            (A - (2/(ntheta+1)) * ((Ag * Ag.') / max(g.' * A * g, 1e-30)));

        xc = min(max(xc, 0), theta_max);
        A = mp_symmetrize_pd(A, theta_max);
    end

    theta_history = theta_history(:, 1:it);
    ellip_trace = ellip_trace(1:it);

    if U >= 3
        thetaSeeds = {};
        thetaSeeds{end+1} = xc(:);
        thetaSeeds{end+1} = 0.75 * xc(:);
        thetaSeeds{end+1} = 1.25 * xc(:);

        for uu = 1:U
            e = zeros(U,1);
            e(uu) = 1;
            thetaSeeds{end+1} = xc(:) + 0.25 * max(1, norm(xc,inf)) * e;
        end

        for uu = 1:U
            e = zeros(U,1);
            e(uu) = 1;
            thetaSeeds{end+1} = max(0, xc(:) - 0.25 * max(1, norm(xc,inf)) * e);
        end

        for ss = 1:numel(thetaSeeds)
            th = min(max(thetaSeeds{ss}, 0), theta_max);
            orderSet_extra = mp_generate_order_candidates(th, atoms, U, opts, prev_best);

            for kk = 1:numel(orderSet_extra)
                cand = mp_inner_solve_for_orders(G, bu_min, w, cb, opts, atoms, orderSet_extra{kk}, th, prev_best);
                [vertex_pool, ~] = mp_pool_add_vertex(vertex_pool, cand, opts.poolUniqueTol);
            end
        end
    end

    used_face_mode = false;
    face_pool = mp_init_vertex_pool(K*N, U);

    if opts.useU2ThetaSweep && U == 2
        rho_grid = logspace(log10(opts.u2ThetaSweepMin), ...
                            log10(opts.u2ThetaSweepMax), ...
                            opts.u2ThetaSweepN);

        if opts.verbose
            fprintf('\nU=2 theta-ratio sweep: %d rho values\n', numel(rho_grid));
        end

        for ii = 1:numel(rho_grid)
            th = [1; rho_grid(ii)];
            orderSet_sweep = mp_generate_order_candidates(th, atoms, U, opts, []);

            for kk = 1:numel(orderSet_sweep)
                cand = mp_inner_solve_for_orders(G, bu_min, w, cb, opts, ...
                    atoms, orderSet_sweep{kk}, th, []);
                [vertex_pool, ~] = mp_pool_add_vertex(vertex_pool, cand, opts.poolUniqueTol);

            end
        end
    end

    if opts.useU2TargetedFinal && U == 2
        orderSet_target = mp_generate_u2_targeted_order_bank(atoms, opts);

        if opts.verbose
            fprintf('\nU=2 joint targeted final harvest: %d order sets\n', numel(orderSet_target));
        end

        for kk = 1:numel(orderSet_target)
            cand_t = mp_inner_solve_targeted(G, bu_min, bu_min, w, cb, opts, ...
                atoms, orderSet_target{kk});
            [vertex_pool, ~] = mp_pool_add_vertex(vertex_pool, cand_t, opts.poolUniqueTol);

            if opts.verbose
                fprintf('  targeted %2d/%2d: E=%-.6g, bu=[%s], short=[%s]\n', ...
                    kk, numel(orderSet_target), cand_t.E, ...
                    num2str(cand_t.bu(:).', ' %.5g'), ...
                    num2str((bu_min(:)-cand_t.bu(:)).', ' %.3e'));
            end
        end
    end

    if opts.useU2ReceiverMACHarvest && U == 2
        mac_pool = mp_harvest_u2_receiver_mac_vertices(G, bu_min, w, cb, opts, atoms);
        if opts.verbose
            fprintf('\nU=2 receiver-MAC diagnostic harvest: %d points\n', size(mac_pool.bu,1));
            for ii = 1:size(mac_pool.bu,1)
                fprintf('  rxMAC %2d: E=%-.6g, bu=[%s]\n', ii, mac_pool.E(ii), ...
                    num2str(mac_pool.bu(ii,:), ' %.5g'));
            end
        end
        for ii = 1:size(mac_pool.bu,1)
            ev = struct();
            ev.bu = mac_pool.bu(ii,:).';
            ev.bu_support = mac_pool.bu_support(ii,:).';
            ev.E = mac_pool.E(ii);
            ev.E_users = mac_pool.E_users(ii,:).';
            ev.powers = reshape(mac_pool.powers(ii,:).', K, N);
            ev.atom_rates = reshape(mac_pool.atom_rates(ii,:).', K, N);
            ev.atom_rates_support = reshape(mac_pool.atom_rates_support(ii,:).', K, N);
            ev.orders = mac_pool.meta{ii}.orders;
            ev.active_rx = mac_pool.meta{ii}.active_rx;
            ev.rx_rate = mac_pool.meta{ii}.rx_rate;
            [vertex_pool, ~] = mp_pool_add_vertex(vertex_pool, ev, opts.poolUniqueTol);
        end
    end

    if opts.enableFaceMode && mp_has_tied_theta_cluster(xc, opts.tieTol)
        face_pool = mp_harvest_tied_theta_face(G, bu_min, w, cb, opts, atoms, xc, prev_best);
        face_pool = mp_pool_lower_envelope(face_pool, opts);
        
        [sel_face, ok_face] = mp_global_select(face_pool, bu_min, opts.hitTol, opts.degenerateTol, opts);
        if ok_face > 0
            sel = sel_face;
            FEAS_FLAG = ok_face;
            if FEAS_FLAG == 1
                sel.mode = 'single-vertex-face';
            else
                sel.mode = 'time-share-face';
            end
            used_face_mode = true;
        end
    end

    raw_vertex_pool = vertex_pool;

    if ~used_face_mode
        [sel, FEAS_FLAG] = mp_global_select(raw_vertex_pool, bu_min, opts.hitTol, opts.degenerateTol, opts);
    end

    vertex_pool = raw_vertex_pool;
    bu_a = sel.bu(:).';

    info = struct();
    info.mode = sel.mode;
    info.target_requested = bu_min(:);
    info.target_supported = mp_tern(FEAS_FLAG > 0, bu_min(:), []);
    info.bu_a = bu_a(:);

    info.E = sel.E;
    info.E_users = sel.E_users(:);
    info.powers = reshape(sel.powers(:), K, N);
    info.atom_energies = info.powers;
    info.atom_rates = reshape(sel.atom_rates(:), K, N);
    info.subuser_rates = info.atom_rates;
    info.atom_rates_support = reshape(sel.atom_rates_support(:), K, N);
    info.atom_table = atoms;

    info.theta_max = theta_max(:);
    info.theta0 = theta0(:);
    info.theta_final = xc(:);
    info.theta_history = theta_history;
    info.ellip_trace = ellip_trace;

    info.vertex_pool_bu = vertex_pool.bu;
    info.vertex_pool_bu_support = vertex_pool.bu_support;
    info.vertex_pool_E = vertex_pool.E(:);
    info.vertex_pool_meta = {vertex_pool.meta{:}};
    pruned_pool = mp_pool_lower_envelope(vertex_pool, opts);
    info.vertex_pool_lower_bu = pruned_pool.bu;
    info.vertex_pool_lower_E = pruned_pool.E(:);
    info.vertex_pool_lower_meta = {pruned_pool.meta{:}};

    info.face_pool_bu = face_pool.bu;
    info.face_pool_bu_support = face_pool.bu_support;
    info.face_pool_E = face_pool.E(:);
    info.face_pool_meta = {face_pool.meta{:}};

    info.best_vertex = sel.best_vertex;
    info.ts_exact = sel.ts_exact;
    info.ts_feasible = sel.ts_feasible;
    info.lp_lambda = sel.lp_lambda;
    info.lp_active_idx = sel.lp_active_idx;
    info.lp_exitflag = sel.lp_exitflag;
    info.lp_message  = sel.lp_message;
    info.G = G;
    info.H = H;
    info.U = U;
    info.K = K;
    info.N = N;
    info.used_face_mode = used_face_mode;

    info.schedule = mp_build_schedule_from_selection(sel, vertex_pool, face_pool, used_face_mode, opts);

    if opts.verbose
        fprintf('\nFinal result:\n');
        fprintf('  FEAS_FLAG = %d\n', FEAS_FLAG);
        fprintf('  bu_a      = [%s]\n', num2str(bu_a, ' %.6g'));
        fprintf('  mode      = %s\n', info.mode);
        fprintf('  E_total   = %.6g\n', info.E);
        fprintf('  E_users   = [%s]\n', num2str(info.E_users.', ' %.6g'));
        fprintf('  theta     = [%s]\n', num2str(info.theta_final.', ' %.6g'));
    end
end

function opts = mp_set_defaults(opts)
    if ~isfield(opts,'useU2ReceiverMACHarvest'), opts.useU2ReceiverMACHarvest = false; end
    if ~isfield(opts,'u2ReceiverMACMaxStarts'), opts.u2ReceiverMACMaxStarts = 12; end
    if ~isfield(opts,'u2ReceiverMACPenalty'), opts.u2ReceiverMACPenalty = 1e8; end
    if ~isfield(opts,'useU2ThetaSweep'), opts.useU2ThetaSweep = false; end
    if ~isfield(opts,'u2ThetaSweepN'),    opts.u2ThetaSweepN = 401; end
    if ~isfield(opts,'u2ThetaSweepMin'),  opts.u2ThetaSweepMin = 1e-3; end
    if ~isfield(opts,'u2ThetaSweepMax'),  opts.u2ThetaSweepMax = 1e3; end
    if ~isfield(opts,'useU2TargetedFinal'), opts.useU2TargetedFinal = false; end
    if ~isfield(opts,'targetPenalty'), opts.targetPenalty = 1e8; end
    if ~isfield(opts,'u2TargetedUseAllLocalOrders'), opts.u2TargetedUseAllLocalOrders = true; end
    if ~isfield(opts,'u2TargetedMaxOrders'), opts.u2TargetedMaxOrders = 36; end
    if ~isfield(opts,'lambdaReportTol'), opts.lambdaReportTol = 1e-2; end
    if ~isfield(opts, 'verbose'), opts.verbose = 1; end
    if ~isfield(opts, 'maxEllipIter'), opts.maxEllipIter = 80; end
    if ~isfield(opts, 'hitTol'), opts.hitTol = 1e-6; end
    if ~isfield(opts, 'tieTol'), opts.tieTol = 5e-3; end
    if ~isfield(opts, 'poolUniqueTol'), opts.poolUniqueTol = 1e-8; end
    if ~isfield(opts, 'degenerateTol'), opts.degenerateTol = 1e-7; end
    if ~isfield(opts, 'nInnerStarts'), opts.nInnerStarts = 8; end
    if ~isfield(opts, 'innerMaxIter'), opts.innerMaxIter = 3000; end
    if ~isfield(opts, 'innerMaxFunEvals'), opts.innerMaxFunEvals = 12000; end
    if ~isfield(opts, 'innerTolX'), opts.innerTolX = 1e-8; end
    if ~isfield(opts, 'innerTolFun'), opts.innerTolFun = 1e-8; end
    if ~isfield(opts, 'bumpBits'), opts.bumpBits = 1; end
    if ~isfield(opts, 'thetaMaxOverride'), opts.thetaMaxOverride = []; end
    if ~isfield(opts, 'thetaInflation'), opts.thetaInflation = 10.0; end
    if ~isfield(opts, 'thetaFloor'), opts.thetaFloor = 1e-8; end
    if ~isfield(opts, 'thetaAdaptiveExpand'), opts.thetaAdaptiveExpand = true; end
    if ~isfield(opts, 'thetaExpandFactor'), opts.thetaExpandFactor = 2.0; end
    if ~isfield(opts, 'thetaNearMaxFrac'), opts.thetaNearMaxFrac = 0.98; end
    if ~isfield(opts, 'thetaExpandShortTol'), opts.thetaExpandShortTol = 1e-3; end
    if ~isfield(opts, 'thetaMaxCap'), opts.thetaMaxCap = Inf; end    % do not cap by default; cap must never shrink theta_max
    if ~isfield(opts, 'capAtTarget'), opts.capAtTarget = true; end
    if ~isfield(opts, 'enableFaceMode'), opts.enableFaceMode = true; end
    if ~isfield(opts, 'faceUniqueTol'), opts.faceUniqueTol = 1e-7; end
    if ~isfield(opts, 'useAvgSupportForUge3'), opts.useAvgSupportForUge3 = false; end
    if ~isfield(opts, 'maxTiePermsPerRx'), opts.maxTiePermsPerRx = 24; end
    if ~isfield(opts, 'maxGlobalOrderCombos'), opts.maxGlobalOrderCombos = 200; end
    if ~isfield(opts, 'nInnerStarts'), opts.nInnerStarts = 20; end
    if ~isfield(opts, 'orderScoreTol'), opts.orderScoreTol = 1e-6; end
    if ~isfield(opts, 'nRandomOrderSweeps'), opts.nRandomOrderSweeps = 12; end
    if ~isfield(opts, 'orderScoreEps'), opts.orderScoreEps = 1e-8; end
    if ~isfield(opts, 'orderPrevMix'), opts.orderPrevMix = 0.5; end
    if ~isfield(opts, 'lpTol'), opts.lpTol = 1e-9; end
    if ~isfield(opts, 'poolDomTol'), opts.poolDomTol = 1e-9; end
    if ~isfield(opts, 'enableHullPrune'), opts.enableHullPrune = true; end
    if ~isfield(opts, 'maxHullPruneVertices'), opts.maxHullPruneVertices = 250; end
    if ~isfield(opts, 'hullPruneEvery'), opts.hullPruneEvery = 25; end
    if ~isfield(opts, 'compressActiveSet'), opts.compressActiveSet = true; end
    if ~isfield(opts, 'useParfor'), opts.useParfor = false; end
    if ~isfield(opts, 'perToneTargetForStarts'), opts.perToneTargetForStarts = true; end
end

function atoms = mp_build_atom_table(U)
% IC atomization: one atom (u,S) for each user u and each nonempty
% receiver subset S that includes u.  Thus user u's information is always
% decodable at its own receiver.
%
% Total K = U * 2^(U-1).
%
% For U=2 this gives:
%   u1_priv, u1_to_12,
%   u2_priv, u2_to_12
%
% Cross-only atoms such as u1_to_2 or u2_to_1 are NOT IC atoms here.

    atoms = repmat(struct( ...
        'index', [], ...
        'label', '', ...
        'user', [], ...
        'decoded_at', [], ...
        'home_rx', [], ...
        'mate_rx', [], ...
        'type', ''), U * 2^(U-1), 1);

    k = 0;

    for u = 1:U
        for mask = 1:(2^U - 1)

            decset = find(bitget(mask, 1:U));
            decset = decset(:).';

            % IC restriction: own receiver must be included.
            if ~ismember(u, decset)
                continue;
            end

            k = k + 1;
            atoms(k).index = k;
            atoms(k).user = u;
            atoms(k).home_rx = u;
            atoms(k).decoded_at = decset;
            atoms(k).mate_rx = setdiff(decset, u);

            if isequal(decset, u)
                atoms(k).type = 'priv';
                atoms(k).label = sprintf('u%d_priv', u);
            else
                atoms(k).type = 'subset';
                atoms(k).label = sprintf('u%d_to_%s', ...
                    u, sprintf('%d', decset));
            end
        end
    end

    atoms = atoms(1:k);
end

function idx = mp_user_atom_indices(atoms, u)
    idx = find([atoms.user] == u);
end

function idx = mp_decodable_atoms_at_receiver(atoms, r)
    K = numel(atoms);
    tf = false(K,1);
    for k = 1:K
        tf(k) = any(atoms(k).decoded_at == r);
    end
    idx = find(tf).';
end

function pool = mp_init_vertex_pool(K, U)
    pool = struct();
    pool.bu = zeros(0, U);
    pool.bu_support = zeros(0, U);
    pool.E = zeros(0, 1);
    pool.E_users = zeros(0, U);
    pool.powers = zeros(0, K);
    pool.atom_rates = zeros(0, K);
    pool.atom_rates_support = zeros(0, K);
    pool.meta = {};
    pool.n_add = 0;
end


function [pool, idx] = mp_pool_add_vertex(pool, ev, tolU)
    idx = [];

    if isempty(pool.bu)
        pool = mp_pool_append_vertex(pool, ev);
        idx = 1;
        return;
    end

    % Keep only the lowest-energy realization for essentially identical b.
    d = sqrt(sum((pool.bu - ev.bu(:).').^2, 2));
    same = find(d <= tolU, 1, 'first');

    if ~isempty(same)
        if ev.E + 1e-14 < pool.E(same)
            pool = mp_pool_replace_vertex(pool, same, ev);
        end
        idx = same;
        return;
    end

    % Do not prune the raw harvested pool online. The previous revision
    % became too aggressive for U>2 and removed the very diversity needed
    % for the final time-sharing LP.
    pool = mp_pool_append_vertex(pool, ev);
    idx = size(pool.bu,1);
end

function pool = mp_pool_append_vertex(pool, ev)
    pool.bu(end+1,:) = ev.bu(:).';
    pool.bu_support(end+1,:) = ev.bu_support(:).';
    pool.E(end+1,1) = ev.E;
    pool.E_users(end+1,:) = ev.E_users(:).';
    pool.powers(end+1,:) = ev.powers(:).';
    pool.atom_rates(end+1,:) = ev.atom_rates(:).';
    pool.atom_rates_support(end+1,:) = ev.atom_rates_support(:).';
    pool.meta{end+1,1} = mp_make_meta(ev);
end

function pool = mp_pool_replace_vertex(pool, i, ev)
    pool.bu(i,:) = ev.bu(:).';
    pool.bu_support(i,:) = ev.bu_support(:).';
    pool.E(i) = ev.E;
    pool.E_users(i,:) = ev.E_users(:).';
    pool.powers(i,:) = ev.powers(:).';
    pool.atom_rates(i,:) = ev.atom_rates(:).';
    pool.atom_rates_support(i,:) = ev.atom_rates_support(:).';
    pool.meta{i} = mp_make_meta(ev);
end

function tf = mp_point_dominates(b1, E1, b2, E2, tol)
    better_or_equal_rates = all(b1(:) >= b2(:) - tol);
    no_worse_energy = (E1 <= E2 + tol);
    strictly_better_somewhere = any(b1(:) > b2(:) + tol) || (E1 < E2 - tol);
    tf = better_or_equal_rates && no_worse_energy && strictly_better_somewhere;
end

function pool = mp_pool_filter(pool, keep)
    keep = keep(:);
    pool.bu = pool.bu(keep,:);
    pool.bu_support = pool.bu_support(keep,:);
    pool.E = pool.E(keep,:);
    pool.E_users = pool.E_users(keep,:);
    pool.powers = pool.powers(keep,:);
    pool.atom_rates = pool.atom_rates(keep,:);
    pool.atom_rates_support = pool.atom_rates_support(keep,:);
    pool.meta = pool.meta(keep,:);
end

function pool = mp_pool_lower_envelope(pool, opts)
    if isempty(pool.bu) || size(pool.bu,1) <= 2
        return;
    end
    if ~isfield(opts,'enableHullPrune') || ~opts.enableHullPrune
        return;
    end
    if ~license('test','optimization_toolbox')
        return;
    end

    nV = size(pool.bu,1);

    % This pruning routine solves one LP per vertex.  For large harvested
    % pools it is only an optional speed/cleanup step and can dominate runtime.
    % Skip it for large pools; mp_global_select still solves the final LP on
    % the raw pool if needed.
    if isfield(opts,'maxHullPruneVertices') && nV > opts.maxHullPruneVertices
        return;
    end

    keep = true(nV,1);
    tol = max(getfield(opts,'lpTol'), 1e-9);

    % A point i is discardable if the remaining points can achieve
    % b >= b_i with weighted energy <= E_i.
    for i = 1:nV
        idx = find(keep);
        idx(idx == i) = [];
        if isempty(idx)
            continue;
        end

        B = pool.bu(idx,:);
        Evec = pool.E(idx);
        target_i = pool.bu(i,:).';
        Ei = pool.E(i);

        f = Evec(:);
        Aineq = -B.';
        bineq = -target_i;
        Aeq = ones(1, numel(idx));
        beq = 1;
        lb = zeros(numel(idx),1);
        ub = ones(numel(idx),1);

        try
            lpOpts = optimoptions('linprog', ...
                'Display','none', ...
                'OptimalityTolerance', tol, ...
                'ConstraintTolerance', tol);
            [lam, val, exitflag] = linprog(f, Aineq, bineq, Aeq, beq, lb, ub, lpOpts);
            if exitflag > 0
                bi = B.' * lam;
                if all(bi >= target_i - 10*tol) && val <= Ei + 10*tol
                    keep(i) = false;
                end
            end
        catch
        end
    end

    pool = mp_pool_filter(pool, keep);
end

function val = mp_getfield_local(s, fld, defaultVal)
    if isstruct(s) && isfield(s, fld)
        val = s.(fld);
    else
        val = defaultVal;
    end
end

function meta = mp_make_meta(ev)
    meta = struct();
    meta.orders = ev.orders;
    meta.powers = ev.powers;
    meta.atom_rates = ev.atom_rates;
    meta.atom_rates_support = ev.atom_rates_support;
    meta.bu = ev.bu(:);
    meta.bu_support = ev.bu_support(:);
    meta.E = ev.E;
    meta.E_users = ev.E_users(:);
    meta.active_rx = ev.active_rx;
    meta.rx_rate = ev.rx_rate;
end

%==========================================================================
% THETA BOUNDS / STARTS / INNER SOLVES
%==========================================================================

function theta_max = mp_compute_theta_max_bounds(G, bu_min, w, cb, opts, atoms)
    U = numel(bu_min);
    theta_max = zeros(U,1);

    theta_probe = ones(U,1);
    probeOrders = mp_generate_order_candidates(theta_probe, atoms, U, opts, []);

    for u = 1:U
        bt = bu_min;
        bt(u) = bt(u) + opts.bumpBits;

        bestE = inf;
        for kk = 1:numel(probeOrders)
            cand = mp_inner_solve_targeted(G, bt, bu_min, w, cb, opts, atoms, probeOrders{kk});
            if cand.E < bestE
                bestE = cand.E;
            end
        end

        theta_max(u) = opts.thetaInflation * max(bestE, opts.thetaFloor);
    end
end

function starts = mp_build_inner_starts(G, bu_min, cb, opts, atoms, prev)
    U = numel(bu_min);
    K = numel(atoms);

    starts = zeros(K, opts.nInnerStarts);
    Pseed = zeros(K,1);

    for u = 1:U
        au = 2^(cb * bu_min(u)) - 1;
        p_iso = max(au / max(G(u,u), eps), 1e-10);
        idxu = mp_user_atom_indices(atoms, u);
        Pseed(idxu) = p_iso / numel(idxu);
    end

    starts(:,1) = log(max(Pseed, 1e-12));

    kk = 2;
    if ~isempty(prev) && kk <= opts.nInnerStarts
        starts(:,kk) = log(max(prev.powers(:), 1e-12));
        kk = kk + 1;
    end

    while kk <= opts.nInnerStarts
        scale = 10.^(0.8 * randn(K,1));
        starts(:,kk) = log(max(Pseed .* scale, 1e-12));
        kk = kk + 1;
    end
end


%==========================================================================
% MULTITONE INNER SOLVES
%==========================================================================

function cand = mp_inner_solve_targeted(G, bt, bu_min, w, cb, opts, atoms, orders)
    if ndims(G) < 3 || size(G,3) == 1
        if ndims(G) == 3
            G1 = G(:,:,1);
        else
            G1 = G;
        end
        cand = mp_inner_solve_targeted_1tone(G1, bt, bu_min, w, cb, opts, atoms, orders);
        return;
    end

    % Multitone targeted solve must enforce the AGGREGATE target across tones.
    % The earlier version incorrectly used bt/N on each tone independently,
    % which can miss a low-energy multitone loading that allocates rates unevenly.
    N = size(G,3);
    K = numel(atoms);
    U = numel(bu_min);

    starts = mp_build_multitone_target_starts(G, bt, bu_min, cb, opts, atoms);
    best = [];

    for ss = 1:size(starts,2)
        z0 = starts(:,ss);
        f = @(zz) mp_targeted_objective_mt(zz, G, bt, w, cb, opts, atoms, orders);

        z = fminsearch(f, z0, optimset( ...
            'Display', 'off', ...
            'TolX', opts.innerTolX, ...
            'TolFun', opts.innerTolFun, ...
            'MaxIter', max(opts.innerMaxIter, 5000), ...
            'MaxFunEvals', max(opts.innerMaxFunEvals, 30000)));

        P = reshape(exp(z(:)), K, N);
        ev = mp_evaluate_fixed_orders_mt(P, G, w, cb, opts, atoms, orders);
        miss = max(bt(:) - ev.bu(:), 0);
        score = ev.E + opts.targetPenalty * sum(miss.^2);

        if isempty(best) || score < best.score
            best = ev;
            best.score = score;
        end
    end

    cand = best;
end

function starts = mp_build_multitone_target_starts(G, bt, bu_min, cb, opts, atoms)
    N = size(G,3);
    K = numel(atoms);
    U = numel(bu_min);
    nS = opts.nInnerStarts;
    starts = zeros(K*N, nS);

    % Baseline: distribute each user's target uniformly over tones and over
    % that user's atoms. This is only a seed; the objective is aggregate.
    Pseed = zeros(K,N);
    for n = 1:N
        for u = 1:U
            au = 2^(cb * max(bt(u),0) / N) - 1;
            p_iso = max(au / max(G(u,u,n), eps), 1e-12);
            idxu = mp_user_atom_indices(atoms, u);
            Pseed(idxu,n) = p_iso / numel(idxu);
        end
    end

    starts(:,1) = log(max(Pseed(:), 1e-12));

    % User-isolated seed: put each user's bits on its direct-channel-best tones.
    if nS >= 2
        Pbest = zeros(K,N);
        for u = 1:U
            [~, ord] = sort(squeeze(G(u,u,:)), 'descend');
            nUse = max(1, min(N, ceil(N/2)));
            use = ord(1:nUse);
            au = 2^(cb * max(bt(u),0) / nUse) - 1;
            idxu = mp_user_atom_indices(atoms, u);
            for nn = use(:).'
                p_iso = max(au / max(G(u,u,nn), eps), 1e-12);
                Pbest(idxu,nn) = p_iso / numel(idxu);
            end
        end
        starts(:,2) = log(max(Pbest(:), 1e-12));
    end

    for ss = 3:nS
        scale = 10.^(0.8 * randn(K,N));
        starts(:,ss) = log(max(Pseed(:) .* scale(:), 1e-12));
    end
end

function val = mp_targeted_objective_mt(z, G, bt, w, cb, opts, atoms, orders)
    N = size(G,3);
    K = numel(atoms);
    P = reshape(exp(z(:)), K, N);
    ev = mp_evaluate_fixed_orders_mt(P, G, w, cb, opts, atoms, orders);
    miss = max(bt(:) - ev.bu(:), 0);
    val = ev.E + opts.targetPenalty * sum(miss.^2) + 1e-12 * sum(P(:).^2);
end

function cand = mp_evaluate_fixed_orders_mt(P, G, w, cb, opts, atoms, orders)
    N = size(G,3);
    K = numel(atoms);
    U = size(G,1);
    tone = cell(N,1);
    for n = 1:N
        tone{n} = mp_evaluate_fixed_orders(P(:,n), G(:,:,n), w, cb, opts, atoms, orders);
    end
    cand = mp_aggregate_tone_candidates(tone, K, U, N, orders);
end

function cand = mp_inner_solve_for_orders(G, bu_min, w, cb, opts, atoms, orders, theta, prev)
    if ndims(G) < 3
        cand = mp_inner_solve_for_orders_1tone(G, bu_min, w, cb, opts, atoms, orders, theta, prev);
        return;
    end

    N = size(G,3);
    K = numel(atoms);
    U = numel(bu_min);

    % The separable Lagrangian is sum_n E_n - theta' * sum_n b_n.
    % Do not cap each tone independently; cap only after aggregation.
    optsTone = opts;
    optsTone.capAtTarget = false;

    % Start-building in the one-tone solver should use a per-tone rate scale.
    bu_tone = bu_min(:) / N;

    tone = cell(N,1);

    if isfield(opts,'useParfor') && opts.useParfor && license('test','Distrib_Computing_Toolbox')
        parfor n = 1:N
            prevn = mp_prev_for_tone(prev, n, K, N);
            tone{n} = mp_inner_solve_for_orders_1tone(G(:,:,n), bu_tone, w, cb, optsTone, atoms, orders, theta, prevn);
        end
    else
        for n = 1:N
            prevn = mp_prev_for_tone(prev, n, K, N);
            tone{n} = mp_inner_solve_for_orders_1tone(G(:,:,n), bu_tone, w, cb, optsTone, atoms, orders, theta, prevn);
        end
    end

    cand = mp_aggregate_tone_candidates(tone, K, U, N, orders);

    if opts.capAtTarget
        reward = min(cand.bu, bu_min(:));
    else
        reward = cand.bu;
    end
    cand.lag_obj = cand.E - theta(:).' * reward(:);
end

function prevn = mp_prev_for_tone(prev, n, K, N)
    prevn = [];
    if isempty(prev) || ~isfield(prev,'powers') || isempty(prev.powers)
        return;
    end
    P = prev.powers;
    if isvector(P)
        if numel(P) == K*N
            P = reshape(P(:), K, N);
        elseif numel(P) == K
            P = P(:);
        else
            return;
        end
    end
    if size(P,2) >= n
        prevn = prev;
        prevn.powers = P(:,n);
    end
end

function cand = mp_aggregate_tone_candidates(tone, K, U, N, orders)
    P = zeros(K,N);
    R = zeros(K,N);
    Rsup = zeros(K,N);
    rx_rate = zeros(K,U,N);
    active_rx = cell(K,N);
    bu = zeros(U,1);
    bu_support = zeros(U,1);
    E_users = zeros(U,1);
    E = 0;

    for n = 1:N
        ev = tone{n};
        P(:,n) = ev.powers(:);
        R(:,n) = ev.atom_rates(:);
        Rsup(:,n) = ev.atom_rates_support(:);
        rx_rate(:,:,n) = ev.rx_rate;
        active_rx(:,n) = ev.active_rx(:);
        bu = bu + ev.bu(:);
        bu_support = bu_support + ev.bu_support(:);
        E_users = E_users + ev.E_users(:);
        E = E + ev.E;
    end

    cand = struct();
    cand.bu = bu(:);
    cand.bu_support = bu_support(:);
    cand.atom_rates = R;
    cand.atom_rates_support = Rsup;
    cand.E = E;
    cand.E_users = E_users(:);
    cand.powers = P;
    cand.orders = orders;
    cand.active_rx = active_rx;
    cand.rx_rate = rx_rate;
end

function cand = mp_inner_solve_targeted_1tone(G, bt, bu_min, w, cb, opts, atoms, orders)
    starts = mp_build_inner_starts(G, bu_min, cb, opts, atoms, []);
    best = [];

    for ss = 1:size(starts,2)
        z0 = starts(:,ss);
        f = @(zz) mp_targeted_objective(zz, G, bt, w, cb, opts, atoms, orders);

        z = fminsearch(f, z0, optimset( ...
            'Display', 'off', ...
            'TolX', opts.innerTolX, ...
            'TolFun', opts.innerTolFun, ...
            'MaxIter', opts.innerMaxIter, ...
            'MaxFunEvals', opts.innerMaxFunEvals));

        P = exp(z(:));
        ev = mp_evaluate_fixed_orders(P, G, w, cb, opts, atoms, orders);
        miss = max(bt(:) - ev.bu(:), 0);
        score = ev.E + 1e6 * sum(miss.^2);

        if isempty(best) || score < best.score
            best = ev;
            best.score = score;
        end
    end

    cand = best;
end

function val = mp_targeted_objective(z, G, bt, w, cb, opts, atoms, orders)
    P = exp(z(:));
    ev = mp_evaluate_fixed_orders(P, G, w, cb, opts, atoms, orders);
    miss = max(bt(:) - ev.bu(:), 0);
    val = ev.E + 1e6 * sum(miss.^2) + 1e-12 * sum(P.^2);
end

function cand = mp_inner_solve_for_orders_1tone(G, bu_min, w, cb, opts, atoms, orders, theta, prev)
    starts = mp_build_inner_starts(G, bu_min, cb, opts, atoms, prev);
    best = [];

    for ss = 1:size(starts,2)
        z0 = starts(:,ss);
        f = @(zz) mp_lagrangian_objective(zz, theta, G, bu_min, w, cb, opts, atoms, orders);

        z = fminsearch(f, z0, optimset( ...
            'Display', 'off', ...
            'TolX', opts.innerTolX, ...
            'TolFun', opts.innerTolFun, ...
            'MaxIter', opts.innerMaxIter, ...
            'MaxFunEvals', opts.innerMaxFunEvals));

        P = exp(z(:));
        ev = mp_evaluate_fixed_orders(P, G, w, cb, opts, atoms, orders);

        if opts.capAtTarget
            reward = min(ev.bu, bu_min);
        else
            reward = ev.bu;
        end

        ev.lag_obj = ev.E - theta(:).' * reward(:);

        if isempty(best) || ev.lag_obj < best.lag_obj
            best = ev;
        end
    end

    cand = best;
end

function val = mp_lagrangian_objective(z, theta, G, bu_min, w, cb, opts, atoms, orders)
    P = exp(z(:));
    ev = mp_evaluate_fixed_orders(P, G, w, cb, opts, atoms, orders);

    if opts.capAtTarget
        reward = min(ev.bu, bu_min);
    else
        reward = ev.bu;
    end

    val = ev.E - theta(:).' * reward(:) + 1e-12 * sum(P.^2);
end

%==========================================================================
% RATE EVALUATION
%==========================================================================

function ev = mp_evaluate_fixed_orders(P, G, w, cb, opts, atoms, orders)
    K = numel(atoms);
    U = size(G,1);
    Cfun = @(x) (1/cb) * log2(1 + max(x,0));

    rx_rate = zeros(K, U);

    for r = 1:U
        ord = orders.rx_orders{r};
        done = false(1, K);
        decmask = false(1, K);
        decmask(ord) = true;

        for m = 1:numel(ord)
            k = ord(m);
            futureDec = decmask & ~done;
            futureDec(k) = false;

            denom = 1;
            for j = 1:K
                uj = atoms(j).user;
                if (~decmask(j)) || futureDec(j)
                    denom = denom + G(r, uj) * P(j);
                end
            end

            num = G(r, atoms(k).user) * P(k);
            rx_rate(k, r) = Cfun(num / denom);
            done(k) = true;
        end
    end

    atom_rates = zeros(K,1);
    atom_rates_support = zeros(K,1);
    bu = zeros(U,1);
    bu_support = zeros(U,1);
    active_rx = cell(K,1);

    for k = 1:K
        decset = atoms(k).decoded_at;
        rr = rx_rate(k, decset).';

        % Exact reliable rate
        rmin = min(rr);
        atom_rates(k) = rmin;

        % Support reward is forced equal to exact reliable rate
        atom_rates_support(k) = rmin;

        % Record active minimizing receivers
        act = decset(abs(rr - rmin) <= 1e-10);
        if isempty(act)
            act = decset(1);
        end
        active_rx{k} = act(:).';

        u = atoms(k).user;
        bu(u) = bu(u) + atom_rates(k);
        bu_support(u) = bu_support(u) + atom_rates_support(k);
    end

    E_users = zeros(U,1);
    for u = 1:U
        idxu = mp_user_atom_indices(atoms, u);
        E_users(u) = sum(P(idxu));
    end
    Et = w(:).' * E_users(:);

    ev = struct();
    ev.bu = bu(:);
    ev.bu_support = bu_support(:);
    ev.atom_rates = atom_rates(:);
    ev.atom_rates_support = atom_rates_support(:);
    ev.E = Et;
    ev.E_users = E_users(:);
    ev.powers = P(:);
    ev.orders = orders;
    ev.active_rx = active_rx;
    ev.rx_rate = rx_rate;
end


%==========================================================================
% U=2 RECEIVER-MAC DIAGNOSTIC HARVEST
%==========================================================================

function pool = mp_harvest_u2_receiver_mac_vertices(G, bu_min, w, cb, opts, atoms)
% Diagnostic/repair harvest for the U=2 dominance/MAC vertices.
%
% This keeps the 4 IC information atoms, but adds candidate vertices obtained
% by viewing each receiver as a two-user MAC between its own private atom and
% the other user's common/cancellable atom.  The rates credited here are the
% receiver-MAC rates at that receiver, not the global Imin rates.  This is
% intended as a diagnostic bridge to test whether the missing exam vertex is
% a receiver-MAC vertex that the current Imin atom evaluator is failing to
% harvest.

    U = 2;
    K = numel(atoms);
    N = size(G,3);
    pool = mp_init_vertex_pool(K*N, U);

    % Receiver 1 MAC atoms: u1_priv and u2_to_12.
    % Receiver 2 MAC atoms: u1_to_12 and u2_priv.
    pairs = { [1 4], [2 3] };

    for r = 1:2
        pair = pairs{r};
        localOrders = [pair; fliplr(pair)];
        for oo = 1:size(localOrders,1)
            ord = localOrders(oo,:);
            cand = mp_inner_solve_u2_receiver_mac(G, bu_min, w, cb, opts, atoms, r, ord);
            [pool, ~] = mp_pool_add_vertex(pool, cand, opts.poolUniqueTol);
        end
    end
end

function cand = mp_inner_solve_u2_receiver_mac(G, bu_min, w, cb, opts, atoms, r, ord)
    K = numel(atoms);
    N = size(G,3);
    nvar = numel(ord)*N;

    starts = mp_build_u2_receiver_mac_starts(G, bu_min, cb, opts, atoms, r, ord);
    nStarts = min(size(starts,2), opts.u2ReceiverMACMaxStarts);
    best = [];

    for ss = 1:nStarts
        z0 = starts(:,ss);
        f = @(zz) mp_u2_receiver_mac_objective(zz, G, bu_min, w, cb, opts, atoms, r, ord);
        z = fminsearch(f, z0, optimset( ...
            'Display','off', ...
            'TolX', opts.innerTolX, ...
            'TolFun', opts.innerTolFun, ...
            'MaxIter', max(opts.innerMaxIter, 5000), ...
            'MaxFunEvals', max(opts.innerMaxFunEvals, 30000)));

        cand = mp_eval_u2_receiver_mac(exp(z(:)), G, w, cb, atoms, r, ord);
        miss = max(bu_min(:) - cand.bu(:), 0);
        score = cand.E + opts.u2ReceiverMACPenalty * sum(miss.^2);
        if isempty(best) || score < best.score
            best = cand;
            best.score = score;
        end
    end

    cand = best;
end

function starts = mp_build_u2_receiver_mac_starts(G, bu_min, cb, opts, atoms, r, ord)
    N = size(G,3);
    nvar = numel(ord)*N;
    nS = max(4, opts.u2ReceiverMACMaxStarts);
    starts = zeros(nvar, nS);

    Pseed = zeros(numel(ord), N);
    for ii = 1:numel(ord)
        k = ord(ii);
        u = atoms(k).user;
        au = 2^(cb * max(bu_min(u),0) / N) - 1;
        for n = 1:N
            Pseed(ii,n) = max(au / max(G(r,u,n), eps), 1e-12);
        end
    end
    starts(:,1) = log(max(Pseed(:), 1e-12));

    % Flat energy seed near the exam's 8/9 per tone scale when N=8.
    Pflat = zeros(numel(ord), N);
    for ii = 1:numel(ord)
        Pflat(ii,:) = max(1/N, 1e-12);
    end
    starts(:,2) = log(max(Pflat(:), 1e-12));

    % User-rate dominant seed: more power to first decoded stream.
    Pdom = Pflat;
    Pdom(1,:) = 2*Pdom(1,:);
    starts(:,3) = log(max(Pdom(:), 1e-12));

    for ss = 4:nS
        scale = 10.^(0.8*randn(size(Pseed)));
        starts(:,ss) = log(max(Pseed(:).*scale(:), 1e-12));
    end
end

function val = mp_u2_receiver_mac_objective(z, G, bu_min, w, cb, opts, atoms, r, ord)
    cand = mp_eval_u2_receiver_mac(exp(z(:)), G, w, cb, atoms, r, ord);
    miss = max(bu_min(:) - cand.bu(:), 0);
    val = cand.E + opts.u2ReceiverMACPenalty * sum(miss.^2) + 1e-12*sum(exp(z(:)).^2);
end

function cand = mp_eval_u2_receiver_mac(pvec, G, w, cb, atoms, r, ord)
    K = numel(atoms);
    U = size(G,1);
    N = size(G,3);
    Cfun = @(x) (1/cb)*log2(1 + max(x,0));

    P = zeros(K,N);
    pmat = reshape(pvec(:), numel(ord), N);
    for ii = 1:numel(ord)
        P(ord(ii),:) = pmat(ii,:);
    end

    R = zeros(K,N);
    Rsup = zeros(K,N);
    rx_rate = zeros(K,U,N);
    active_rx = cell(K,N);
    bu = zeros(U,1);

    for n = 1:N
        done = false(1,K);
        decmask = false(1,K);
        decmask(ord) = true;
        for m = 1:numel(ord)
            k = ord(m);
            futureDec = decmask & ~done;
            futureDec(k) = false;
            denom = 1;
            for jj = 1:numel(ord)
                j = ord(jj);
                if futureDec(j)
                    uj = atoms(j).user;
                    denom = denom + G(r,uj,n)*P(j,n);
                end
            end
            u = atoms(k).user;
            rate = Cfun(G(r,u,n)*P(k,n)/denom);
            rx_rate(k,r,n) = rate;
            R(k,n) = rate;
            Rsup(k,n) = rate;
            active_rx{k,n} = r;
            bu(u) = bu(u) + rate;
            done(k) = true;
        end
    end

    E_users = zeros(U,1);
    for u = 1:U
        idxu = mp_user_atom_indices(atoms,u);
        E_users(u) = sum(P(idxu,:),'all');
    end

    orders = struct('rx_orders',{cell(U,1)});
    orders.rx_orders{r} = ord;
    for rr = setdiff(1:U,r)
        orders.rx_orders{rr} = [];
    end

    cand = struct();
    cand.bu = bu(:);
    cand.bu_support = bu(:);
    cand.atom_rates = R;
    cand.atom_rates_support = Rsup;
    cand.E_users = E_users(:);
    cand.E = w(:).'*E_users(:);
    cand.powers = P;
    cand.orders = orders;
    cand.active_rx = active_rx;
    cand.rx_rate = rx_rate;
end

%==========================================================================
% ORDER GENERATION / FACE HARVEST
%==========================================================================


function orderSet = mp_generate_u2_targeted_order_bank(atoms, opts)
% U=2 diagnostic/repair order bank for the final joint targeted solve.
% This is deliberately small: each receiver has three decodable atoms in
% the 4-atom IC model, so all local permutations give at most 6*6 = 36
% receiver-order sets.  This avoids assuming the same-user private/common
% tie order while keeping the main Section-5.6 theta logic unchanged.

    U = 2;
    local = cell(U,1);
    for r = 1:U
        dec = mp_decodable_atoms_at_receiver(atoms, r);
        if opts.u2TargetedUseAllLocalOrders
            local{r} = perms(dec);
        else
            % Fallback: use the deterministic theta rule's two user priorities.
            Oa = mp_generate_order_candidates([2;1], atoms, U, opts, []);
            Ob = mp_generate_order_candidates([1;2], atoms, U, opts, []);
            orderSet = {Oa{1}, Ob{1}};
            return;
        end
    end

    n1 = size(local{1},1);
    n2 = size(local{2},1);
    nTot = min(n1*n2, opts.u2TargetedMaxOrders);
    orderSet = cell(nTot,1);
    t = 0;
    for i = 1:n1
        for j = 1:n2
            t = t + 1;
            if t > nTot
                return;
            end
            tmp = struct('rx_orders', {cell(U,1)});
            tmp.rx_orders{1} = local{1}(i,:);
            tmp.rx_orders{2} = local{2}(j,:);
            orderSet{t} = tmp;
        end
    end
end

function tf = mp_has_tied_theta_cluster(theta, tieTol)
    tf = false;
    U = numel(theta);
    for i = 1:U-1
        for j = i+1:U
            if abs(theta(i)-theta(j)) <= tieTol
                tf = true;
                return;
            end
        end
    end
end

function orderSet = mp_generate_order_candidates(theta, atoms, U, opts, prev_best) %#ok<INUSD>
% Section-5.6 deterministic theta-to-order rule.
%
% A theta vector determines the user order by relative theta size.
% Larger theta user is decoded earlier.  This produces ONE global
% receiver-order set, except near exact theta ties, where tied user
% permutations may be generated for face/time-sharing diagnostics.

    theta = theta(:);

    % User priority: largest theta first.
    [thetaSorted, userOrder] = sort(theta, 'descend');

    % Detect only true/near theta ties.
    hasTie = false;
    for i = 1:U-1
        if abs(thetaSorted(i) - thetaSorted(i+1)) <= opts.tieTol
            hasTie = true;
            break;
        end
    end

    if ~hasTie
        userOrders = userOrder(:).';
    else
        % For small U this is fine; for larger U cap by maxTiePermsPerRx.
        P = perms(1:U);
        % Keep permutations that only reorder tied users approximately.
        % Simpler safe version: try first few full user permutations.
        nKeep = min(size(P,1), max(1, opts.maxTiePermsPerRx));
        userOrders = P(1:nKeep,:);
    end

    orderSet = cell(size(userOrders,1),1);

    for q = 1:size(userOrders,1)
        uord = userOrders(q,:);

        tmp = struct('rx_orders', {cell(U,1)});

        for r = 1:U
            rxord = [];

            for uu = uord
                dec = mp_decodable_atoms_at_receiver(atoms, r);
                dec_u = dec(arrayfun(@(k) atoms(k).user == uu, dec));

                if isempty(dec_u)
                    continue;
                end

                % Same-user atoms have the same theta.  Use a fixed
                % deterministic convention: more-common atoms first,
                % private/less-common atoms later.
                dec_sizes = arrayfun(@(k) numel(atoms(k).decoded_at), dec_u);
                [~, jj] = sortrows([-dec_sizes(:), dec_u(:)]);
                dec_u = dec_u(jj);

                rxord = [rxord, dec_u(:).']; %#ok<AGROW>
            end

            tmp.rx_orders{r} = rxord;
        end

        orderSet{q} = tmp;
    end
end

function O = mp_local_orders_from_scores(atomIdx, scores, tieTol, maxPerms)
    [ws, sidx] = sort(scores(:), 'ascend');
    atomSorted = atomIdx(sidx(:)).';

    groups = cell(0,1);
    gstart = 1;
    for i = 2:numel(ws)
        if abs(ws(i) - ws(i-1)) > tieTol
            groups{end+1,1} = atomSorted(gstart:i-1); %#ok<AGROW>
            gstart = i;
        end
    end
    groups{end+1,1} = atomSorted(gstart:end);

    blockPerms = cell(numel(groups),1);
    for g = 1:numel(groups)
        grp = groups{g};
        if numel(grp) == 1
            blockPerms{g} = grp;
        else
            Pg = perms(grp);
            if size(Pg,1) > maxPerms
                Pg = Pg(1:maxPerms, :);
            end
            blockPerms{g} = Pg;
        end
    end

    O = blockPerms{1};
    for g = 2:numel(blockPerms)
        A1 = O;
        A2 = blockPerms{g};
        newO = zeros(size(A1,1) * size(A2,1), size(A1,2) + size(A2,2));
        tt = 0;
        for i = 1:size(A1,1)
            for j = 1:size(A2,1)
                tt = tt + 1;
                newO(tt,:) = [A1(i,:), A2(j,:)];
            end
        end
        O = newO;
    end

    O = unique(O, 'rows');
    if size(O,1) > maxPerms
        O = O(1:maxPerms, :);
    end
end

function facePool = mp_harvest_tied_theta_face(G, bu_min, w, cb, opts, atoms, theta, prev_best)
    U = numel(theta);

    theta_face = theta(:);
    orderSet = mp_generate_order_candidates(theta_face, atoms, U, opts, prev_best);

    K = numel(atoms);
    if ndims(G) < 3
        N = 1;
    else
        N = size(G,3);
    end

    % Multitone candidates store K-by-N powers/rates, flattened in the pool.
    % The old code used K here, which crashes for N>1 when ev.powers has K*N entries.
    facePool = mp_init_vertex_pool(K*N, U);

    for kk = 1:numel(orderSet)
        cand = mp_inner_solve_for_orders(G, bu_min, w, cb, opts, atoms, orderSet{kk}, theta_face, prev_best);
        facePool = mp_pool_add_vertex(facePool, cand, opts.faceUniqueTol);
    end
end

%==========================================================================
% FINAL GLOBAL SELECTION
%==========================================================================


function [sel, FEAS_FLAG] = mp_global_select(pool, target, hitTol, degTol, opts)
    rawPool = pool;

    try
        prunedPool = mp_pool_lower_envelope(pool, opts);
        if isempty(prunedPool.bu)
            prunedPool = rawPool;
        end
    catch
        prunedPool = rawPool;
    end

    [sel_p, flag_p] = mp_global_select_core(prunedPool, target, hitTol, degTol, opts);
    [sel_r, flag_r] = mp_global_select_core(rawPool,    target, hitTol, degTol, opts);

    % Prefer feasible over infeasible; if both feasible, prefer lower energy.
    if flag_r > flag_p || ...
       (flag_r > 0 && flag_p > 0 && sel_r.E < sel_p.E - 1e-8)
        sel = sel_r;
        FEAS_FLAG = flag_r;
    else
        sel = sel_p;
        FEAS_FLAG = flag_p;
    end
end

function [sel, FEAS_FLAG] = mp_global_select_core(pool, target, hitTol, degTol, opts)
    target = target(:);
    U = numel(target);
    K = size(pool.powers, 2);
    nV = size(pool.bu, 1);

    FEAS_FLAG = 0;
    sel = struct();
    sel.mode = 'target-miss';
    sel.bu = zeros(U,1);
    sel.E = inf;
    sel.E_users = zeros(U,1);
    sel.powers = zeros(K,1);
    sel.atom_rates = zeros(K,1);
    sel.atom_rates_support = zeros(K,1);
    sel.best_vertex = struct('idx',[],'bu',[],'E',[],'powers',[],'atom_rates',[],'meta',[]);
    sel.ts_exact = struct('used',0,'kind','','frac',[],'idx',[],'bu_vertices',[],'E_vertices',[],'E_final',[]);
    sel.ts_feasible = struct('used',0,'kind','','frac',[],'idx',[],'bu_vertices',[],'E_vertices',[],'E_final',[]);
    sel.lp_lambda = [];
    sel.lp_active_idx = [];
    sel.lp_active_meta = {};
    sel.lp_exitflag = [];
    sel.lp_message = '';

    if nV == 0
        return;
    end

    bestShort = inf;
    bestIdx = 1;
    for i = 1:nV
        vi = pool.bu(i,:).';
        shorti = sum(max(target - vi, 0).^2);
        if shorti < bestShort || (abs(shorti-bestShort)<=1e-14 && pool.E(i) < pool.E(bestIdx))
            bestShort = shorti;
            bestIdx = i;
        end

       
    end

    B = pool.bu;
    Evec = pool.E(:);

    f = Evec;
    Aineq = -B.';
    bineq = -target;
    Aeq = ones(1, nV);
    beq = 1;
    lb = zeros(nV,1);
    ub = ones(nV,1);

    lam = [];
    exitflag = -999;
    lp_message = '';

    try
        if license('test','optimization_toolbox')
            lpOpts = optimoptions('linprog', ...
                'Display', 'none', ...
                'OptimalityTolerance', max(opts.lpTol,1e-9), ...
                'ConstraintTolerance', max(opts.lpTol,1e-9));
            [lam, ~, exitflag, output] = linprog(f, Aineq, bineq, Aeq, beq, lb, ub, lpOpts);
            if isfield(output, 'message')
                lp_message = output.message;
            end
        else
            lam = [];
            exitflag = -999;
            lp_message = 'Optimization Toolbox / linprog not available';
        end
    catch ME
        lam = [];
        exitflag = -999;
        lp_message = ME.message;
    end

    sel.lp_exitflag = exitflag;
    sel.lp_message  = lp_message;

    if ~isempty(lam) && exitflag > 0
        lam = mp_compress_lambda_support(lam, B, Evec, target, opts);
        bu_mix = B.' * lam;
        if all(bu_mix >= target - 10*hitTol)
            FEAS_FLAG = 2;
            sel.mode = 'time-share-lp';
            sel.bu = bu_mix;
            sel.E = Evec.' * lam;
            sel.E_users = (pool.E_users.' * lam);
            sel.powers = (pool.powers.' * lam);
            sel.atom_rates = (pool.atom_rates.' * lam);
            sel.atom_rates_support = (pool.atom_rates_support.' * lam);
            sel.best_vertex = mp_pack_best_vertex(pool, bestIdx);
            sel.lp_lambda = lam;
            sel.lp_active_idx = find(lam > max(opts.lpTol,1e-8));
            sel.lp_active_meta = pool.meta(sel.lp_active_idx);

            sel.ts_feasible.used = 1;
            sel.ts_feasible.kind = 'lp-feasible';
            sel.ts_feasible.frac = lam(:);
            sel.ts_feasible.idx = sel.lp_active_idx(:);
            sel.ts_feasible.bu_vertices = pool.bu(sel.lp_active_idx, :);
            sel.ts_feasible.E_vertices = pool.E(sel.lp_active_idx);
            sel.ts_feasible.E_final = sel.E;
            return;
        end
    end

    sel.mode = 'target-miss';
    sel.bu = pool.bu(bestIdx,:).';
    sel.E = pool.E(bestIdx);
    sel.E_users = pool.E_users(bestIdx,:).';
    sel.powers = pool.powers(bestIdx,:).';
    sel.atom_rates = pool.atom_rates(bestIdx,:).';
    sel.atom_rates_support = pool.atom_rates_support(bestIdx,:).';
    sel.best_vertex = mp_pack_best_vertex(pool, bestIdx);
end

%==========================================================================
% MISC HELPERS%==========================================================================
% MISC HELPERS
%==========================================================================

function A = mp_symmetrize_pd(A, theta_max)
    A = 0.5*(A+A.');
    [V,D] = eig(A);
    d = diag(D);
    d = max(d, 1e-14*max(1,max(theta_max.^2)));
    A = V*diag(d)*V.';
    A = 0.5*(A+A.');
end

function out = mp_pack_best_vertex(pool, idx)
    out = struct();
    out.idx = idx;
    out.bu = pool.bu(idx,:).';
    out.E = pool.E(idx);
    out.powers = pool.powers(idx,:).';
    out.atom_rates = pool.atom_rates(idx,:).';
    out.meta = pool.meta{idx};
end

function ts = mp_pack_ts(pool, i, j, alpha, kind)
    ts = struct();
    ts.used = 1;
    ts.kind = kind;
    ts.frac = alpha(:);
    ts.idx = [i; j];
    ts.bu_vertices = [pool.bu(i,:); pool.bu(j,:)];
    ts.E_vertices = [pool.E(i); pool.E(j)];
    ts.E_final = alpha(1)*pool.E(i) + alpha(2)*pool.E(j);
end

function [tf, alpha] = mp_segment_exact_target(v1, v2, t, tol)
    tf = false;
    alpha = [];

    D = v1 - v2;
    idx = find(abs(D) > tol);
    if isempty(idx)
        if norm(v1 - t, inf) <= tol
            tf = true;
            alpha = [1; 0];
        end
        return;
    end

    vals = zeros(numel(idx),1);
    for k = 1:numel(idx)
        ii = idx(k);
        vals(k) = (t(ii) - v2(ii)) / D(ii);
    end

    a1 = vals(1);
    if max(abs(vals - a1)) > 1e-5
        return;
    end

    a2 = 1 - a1;
    mix = a1*v1 + a2*v2;

    if a1 >= -tol && a2 >= -tol && norm(mix - t, inf) <= 10*tol
        tf = true;
        alpha = [a1; a2];
    end
end

function [tf, lamLo, lamHi] = mp_segment_dominates_target(v1, v2, t, tol)
    U = numel(t);
    lamLo = 0;
    lamHi = 1;

    for k = 1:U
        d = v1(k) - v2(k);
        if abs(d) <= tol
            if v2(k) < t(k) - tol
                tf = false;
                lamLo = NaN; lamHi = NaN;
                return;
            end
        elseif d > 0
            lamLo = max(lamLo, (t(k)-v2(k))/d);
        else
            lamHi = min(lamHi, (t(k)-v2(k))/d);
        end
    end

    lamLo = max(lamLo,0);
    lamHi = min(lamHi,1);
    tf = (lamLo <= lamHi + tol);
end

function lam = mp_closest_lambda_to_target(v1, v2, t, lamLo, lamHi)
    d = v1-v2;
    if norm(d) <= 0
        lam = lamLo;
        return;
    end
    lam0 = (d.'*(t-v2)) / (d.'*d);
    lam = min(max(lam0, lamLo), lamHi);
end


function lam = mp_compress_lambda_support(lam, B, Evec, target, opts)
    lam = lam(:);
    if ~isfield(opts,'compressActiveSet') || ~opts.compressActiveSet
        return;
    end
    if ~license('test','optimization_toolbox')
        return;
    end

    tol = max(opts.lpTol, 1e-9);
    active = find(lam > max(tol,1e-8));
    U = size(B,2);

    % On a proper supporting face in R^U, at most U active points are needed.
    % We greedily test whether one active point can be removed while preserving
    % feasibility and objective value.
    while numel(active) > U
        dropped = false;
        for kk = 1:numel(active)
            idx = active;
            idx(kk) = [];
            Bk = B(idx,:);
            Ek = Evec(idx);

            f = Ek(:);
            Aineq = -Bk.';
            bineq = -target(:);
            Aeq = ones(1, numel(idx));
            beq = 1;
            lb = zeros(numel(idx),1);
            ub = ones(numel(idx),1);

            try
                lpOpts = optimoptions('linprog', ...
                    'Display','none', ...
                    'OptimalityTolerance', tol, ...
                    'ConstraintTolerance', tol);
                [lamk, valk, exitflag] = linprog(f, Aineq, bineq, Aeq, beq, lb, ub, lpOpts);
                if exitflag > 0
                    buk = Bk.' * lamk;
                    if all(buk >= target(:) - 10*tol) && valk <= Evec(:).' * lam + 10*tol
                        lamNew = zeros(size(lam));
                        lamNew(idx) = lamk;
                        lam = lamNew;
                        active = find(lam > max(tol,1e-8));
                        dropped = true;
                        break;
                    end
                end
            catch
            end
        end
        if ~dropped
            break;
        end
    end
end

function y = mp_tern(cond, a, b)
    if cond
        y = a;
    else
        y = b;
    end
end

function sched = mp_build_schedule_from_selection(sel, ~, ~, ~, opts)

    if isfield(sel,'lp_active_meta') && ~isempty(sel.lp_active_meta)
        idxAll  = sel.lp_active_idx(:);
        lamAll  = sel.lp_lambda(idxAll);
        metaAll = sel.lp_active_meta(:);

        keep = lamAll >= opts.lambdaReportTol;

        % If lambdaReportTol removes everything, keep all active shares.
        if ~any(keep)
            keep = true(size(lamAll));
        end

        idx  = idxAll(keep);
        lam  = lamAll(keep);
        meta = metaAll(keep);

        lam = lam / sum(lam);
        nS = numel(lam);
    else
        idx = [];
        lam = 1;
        meta = {};
        nS = 1;
    end

    sched = repmat(struct( ...
        'lambda', [], ...
        'pool_index', [], ...
        'powers', [], ...
        'atom_energies', [], ...
        'atom_rates', [], ...
        'atom_rates_support', [], ...
        'bu', [], ...
        'E', [], ...
        'E_users', [], ...
        'orders', [], ...
        'rx_orders', [], ...
        'rx_rate', [], ...
        'active_rx', []), nS, 1);

    for s = 1:nS
        if isfield(sel,'lp_active_meta') && ~isempty(sel.lp_active_meta)
            m = meta{s};
            sched(s).lambda = lam(s);
            sched(s).pool_index = idx(s);
        else
            m = sel.best_vertex.meta;
            sched(s).lambda = 1;
            sched(s).pool_index = sel.best_vertex.idx;
        end

        sched(s).powers = m.powers;
        sched(s).atom_energies = m.powers;
        sched(s).atom_rates = m.atom_rates;
        sched(s).atom_rates_support = m.atom_rates_support;
        sched(s).bu = m.bu(:);
        sched(s).E = m.E;
        sched(s).E_users = m.E_users(:);
        sched(s).orders = m.orders;
        sched(s).rx_orders = m.orders.rx_orders;
        sched(s).rx_rate = m.rx_rate;
        sched(s).active_rx = m.active_rx;
    end
end